<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>Chapter 18 Cookbook – Cool things to do with it &mdash; Biopython_en 1.0 documentation</title>
    
    <link rel="stylesheet" href="_static/default.css" type="text/css" />
    <link rel="stylesheet" href="_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    './',
        VERSION:     '1.0',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="_static/jquery.js"></script>
    <script type="text/javascript" src="_static/underscore.js"></script>
    <script type="text/javascript" src="_static/doctools.js"></script>
    <link rel="top" title="Biopython_en 1.0 documentation" href="index.html" />
    <link rel="next" title="Chapter 19 The Biopython testing framework" href="chr19.html" />
    <link rel="prev" title="Chapter 17 Graphics including GenomeDiagram" href="chr17.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="chr19.html" title="Chapter 19 The Biopython testing framework"
             accesskey="N">next</a> |</li>
        <li class="right" >
          <a href="chr17.html" title="Chapter 17 Graphics including GenomeDiagram"
             accesskey="P">previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <div class="section" id="chapter-18-cookbook-cool-things-to-do-with-it">
<h1>Chapter 18  Cookbook – Cool things to do with it<a class="headerlink" href="#chapter-18-cookbook-cool-things-to-do-with-it" title="Permalink to this headline">¶</a></h1>
<p>Biopython now has two collections of “cookbook” examples – this chapter
(which has been included in this tutorial for many years and has
gradually grown), and
<tt class="docutils literal"><span class="pre">`http://biopython.org/wiki/Category:Cookbook</span></tt> &lt;<a class="reference external" href="http://biopython.org/wiki/Category:Cookbook">http://biopython.org/wiki/Category:Cookbook</a>&gt;`__
which is a user contributed collection on our wiki.</p>
<p>We’re trying to encourage Biopython users to contribute their own
examples to the wiki. In addition to helping the community, one direct
benefit of sharing an example like this is that you could also get some
feedback on the code from other Biopython users and developers - which
could help you improve all your Python code.</p>
<p>In the long term, we may end up moving all of the examples in this
chapter to the wiki, or elsewhere within the tutorial.</p>
<div class="section" id="working-with-sequence-files">
<h2>18.1  Working with sequence files<a class="headerlink" href="#working-with-sequence-files" title="Permalink to this headline">¶</a></h2>
<p>This section shows some more examples of sequence input/output, using
the <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> module described in
Chapter <a class="reference external" href="#chapter:Bio.SeqIO">5</a>.</p>
<div class="section" id="filtering-a-sequence-file">
<h3>18.1.1  Filtering a sequence file<a class="headerlink" href="#filtering-a-sequence-file" title="Permalink to this headline">¶</a></h3>
<p>Often you’ll have a large file with many sequences in it (e.g. FASTA
file or genes, or a FASTQ or SFF file of reads), a separate shorter list
of the IDs for a subset of sequences of interest, and want to make a new
sequence file for this subset.</p>
<p>Let’s say the list of IDs is in a simple text file, as the first word on
each line. This could be a tabular file where the first column is the
ID. Try something like this:</p>
<p>Note that we use a Python <tt class="docutils literal"><span class="pre">set</span></tt> rather than a <tt class="docutils literal"><span class="pre">list</span></tt>, this makes
testing membership faster.</p>
</div>
<div class="section" id="producing-randomised-genomes">
<h3>18.1.2  Producing randomised genomes<a class="headerlink" href="#producing-randomised-genomes" title="Permalink to this headline">¶</a></h3>
<p>Let’s suppose you are looking at genome sequence, hunting for some
sequence feature – maybe extreme local GC% bias, or possible restriction
digest sites. Once you’ve got your Python code working on the real
genome it may be sensible to try running the same search on randomised
versions of the same genome for statistical analysis (after all, any
“features” you’ve found could just be there just by chance).</p>
<p>For this discussion, we’ll use the GenBank file for the pPCP1 plasmid
from <em>Yersinia pestis biovar Microtus</em>. The file is included with the
Biopython unit tests under the GenBank folder, or you can get it from
our website,
<tt class="docutils literal"><span class="pre">`NC_005816.gb</span></tt> &lt;<a class="reference external" href="http://biopython.org/SRC/biopython/Tests/GenBank/NC_005816.gb">http://biopython.org/SRC/biopython/Tests/GenBank/NC_005816.gb</a>&gt;`__.
This file contains one and only one record, so we can read it in as a
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> using the <tt class="docutils literal"><span class="pre">Bio.SeqIO.read()</span></tt> function:</p>
<p>So, how can we generate a shuffled versions of the original sequence? I
would use the built in Python <tt class="docutils literal"><span class="pre">random</span></tt> module for this, in particular
the function <tt class="docutils literal"><span class="pre">random.shuffle</span></tt> – but this works on a Python list. Our
sequence is a <tt class="docutils literal"><span class="pre">Seq</span></tt> object, so in order to shuffle it we need to turn
it into a list:</p>
<p>Now, in order to use <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to output the shuffled sequence, we
need to construct a new <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> with a new <tt class="docutils literal"><span class="pre">Seq</span></tt> object using
this shuffled list. In order to do this, we need to turn the list of
nucleotides (single letter strings) into a long string – the standard
Python way to do this is with the string object’s join method.</p>
<p>Let’s put all these pieces together to make a complete Python script
which generates a single FASTA file containing 30 randomly shuffled
versions of the original sequence.</p>
<p>This first version just uses a big for loop and writes out the records
one by one (using the <tt class="docutils literal"><span class="pre">SeqRecord</span></tt>’s format method described in
Section <a class="reference external" href="#sec:Bio.SeqIO-and-StringIO">5.5.4</a>):</p>
<p>Personally I prefer the following version using a function to shuffle
the record and a generator expression instead of the for loop:</p>
</div>
<div class="section" id="translating-a-fasta-file-of-cds-entries">
<h3>18.1.3  Translating a FASTA file of CDS entries<a class="headerlink" href="#translating-a-fasta-file-of-cds-entries" title="Permalink to this headline">¶</a></h3>
<p>Suppose you’ve got an input file of CDS entries for some organism, and
you want to generate a new FASTA file containing their protein
sequences. i.e. Take each nucleotide sequence from the original file,
and translate it. Back in Section <a class="reference external" href="#sec:translation">3.9</a> we saw
how to use the <tt class="docutils literal"><span class="pre">Seq</span></tt> object’s <tt class="docutils literal"><span class="pre">translate</span> <span class="pre">method</span></tt>, and the optional
<tt class="docutils literal"><span class="pre">cds</span></tt> argument which enables correct translation of alternative start
codons.</p>
<p>We can combine this with <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> as shown in the reverse
complement example in
Section <a class="reference external" href="#sec:SeqIO-reverse-complement">5.5.3</a>. The key point is
that for each nucleotide <tt class="docutils literal"><span class="pre">SeqRecord</span></tt>, we need to create a protein
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> - and take care of naming it.</p>
<p>You can write you own function to do this, choosing suitable protein
identifiers for your sequences, and the appropriate genetic code. In
this example we just use the default table and add a prefix to the
identifier:</p>
<p>We can then use this function to turn the input nucleotide records into
protein records ready for output. An elegant way and memory efficient
way to do this is with a generator expression:</p>
<p>This should work on any FASTA file of complete coding sequences. If you
are working on partial coding sequences, you may prefer to use
<tt class="docutils literal"><span class="pre">nuc_record.seq.translate(to_stop=True)</span></tt> in the example above, as this
wouldn’t check for a valid start codon etc.</p>
</div>
<div class="section" id="making-the-sequences-in-a-fasta-file-upper-case">
<h3>18.1.4  Making the sequences in a FASTA file upper case<a class="headerlink" href="#making-the-sequences-in-a-fasta-file-upper-case" title="Permalink to this headline">¶</a></h3>
<p>Often you’ll get data from collaborators as FASTA files, and sometimes
the sequences can be in a mixture of upper and lower case. In some cases
this is deliberate (e.g. lower case for poor quality regions), but
usually it is not important. You may want to edit the file to make
everything consistent (e.g. all upper case), and you can do this easily
using the <tt class="docutils literal"><span class="pre">upper()</span></tt> method of the <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object (added in
Biopython 1.55):</p>
<p>How does this work? The first line is just importing the <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt>
module. The second line is the interesting bit – this is a Python
generator expression which gives an upper case version of each record
parsed from the input file (<tt class="docutils literal"><span class="pre">mixed.fas</span></tt>). In the third line we give
this generator expression to the <tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt> function and it
saves the new upper cases records to our output file (<tt class="docutils literal"><span class="pre">upper.fas</span></tt>).</p>
<p>The reason we use a generator expression (rather than a list or list
comprehension) is this means only one record is kept in memory at a
time. This can be really important if you are dealing with large files
with millions of entries.</p>
</div>
<div class="section" id="sorting-a-sequence-file">
<h3>18.1.5  Sorting a sequence file<a class="headerlink" href="#sorting-a-sequence-file" title="Permalink to this headline">¶</a></h3>
<p>Suppose you wanted to sort a sequence file by length (e.g. a set of
contigs from an assembly), and you are working with a file format like
FASTA or FASTQ which <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> can read, write (and index).</p>
<p>If the file is small enough, you can load it all into memory at once as
a list of <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects, sort the list, and save it:</p>
<p>The only clever bit is specifying a comparison function for how to sort
the records (here we sort them by length). If you wanted the longest
records first, you could flip the comparison or use the reverse
argument:</p>
<p>Now that’s pretty straight forward - but what happens if you have a very
large file and you can’t load it all into memory like this? For example,
you might have some next-generation sequencing reads to sort by length.
This can be solved using the <tt class="docutils literal"><span class="pre">Bio.SeqIO.index()</span></tt> function.</p>
<p>First we scan through the file once using <tt class="docutils literal"><span class="pre">Bio.SeqIO.parse()</span></tt>,
recording the record identifiers and their lengths in a list of tuples.
We then sort this list to get them in length order, and discard the
lengths. Using this sorted list of identifiers <tt class="docutils literal"><span class="pre">Bio.SeqIO.index()</span></tt>
allows us to retrieve the records one by one, and we pass them to
<tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt> for output.</p>
<p>These examples all use <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to parse the records into
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects which are output using <tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt>. What
if you want to sort a file format which <tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt> doesn’t
support, like the plain text SwissProt format? Here is an alternative
solution using the <tt class="docutils literal"><span class="pre">get_raw()</span></tt> method added to <tt class="docutils literal"><span class="pre">Bio.SeqIO.index()</span></tt>
in Biopython 1.54 (see Section <a class="reference external" href="#sec:seqio-index-getraw">5.4.2.2</a>).</p>
<p>As a bonus, because it doesn’t parse the data into <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects
a second time it should be faster.</p>
</div>
<div class="section" id="simple-quality-filtering-for-fastq-files">
<h3>18.1.6  Simple quality filtering for FASTQ files<a class="headerlink" href="#simple-quality-filtering-for-fastq-files" title="Permalink to this headline">¶</a></h3>
<p>The FASTQ file format was introduced at Sanger and is now widely used
for holding nucleotide sequencing reads together with their quality
scores. FASTQ files (and the related QUAL files) are an excellent
example of per-letter-annotation, because for each nucleotide in the
sequence there is an associated quality score. Any per-letter-annotation
is held in a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> in the <tt class="docutils literal"><span class="pre">letter_annotations</span></tt> dictionary as a
list, tuple or string (with the same number of elements as the sequence
length).</p>
<p>One common task is taking a large set of sequencing reads and filtering
them (or cropping them) based on their quality scores. The following
example is very simplistic, but should illustrate the basics of working
with quality data in a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object. All we are going to do here
is read in a file of FASTQ data, and filter it to pick out only those
records whose PHRED quality scores are all above some threshold (here
20).</p>
<p>For this example we’ll use some real data downloaded from the ENA
sequence read archive,
<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</a>&gt;`__
(2MB) which unzips to a 19MB file <tt class="docutils literal"><span class="pre">SRR020192.fastq</span></tt>. This is some
Roche 454 GS FLX single end data from virus infected California sea
lions (see
<tt class="docutils literal"><span class="pre">`http://www.ebi.ac.uk/ena/data/view/SRS004476</span></tt> &lt;<a class="reference external" href="http://www.ebi.ac.uk/ena/data/view/SRS004476">http://www.ebi.ac.uk/ena/data/view/SRS004476</a>&gt;`__
for details).</p>
<p>First, let’s count the reads:</p>
<p>Now let’s do a simple filtering for a minimum PHRED quality of 20:</p>
<p>This pulled out only 14580 reads out of the 41892 present. A more
sensible thing to do would be to quality trim the reads, but this is
intended as an example only.</p>
<p>FASTQ files can contain millions of entries, so it is best to avoid
loading them all into memory at once. This example uses a generator
expression, which means only one <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> is created at a time -
avoiding any memory limitations.</p>
</div>
<div class="section" id="trimming-off-primer-sequences">
<h3>18.1.7  Trimming off primer sequences<a class="headerlink" href="#trimming-off-primer-sequences" title="Permalink to this headline">¶</a></h3>
<p>For this example we’re going to pretend that <tt class="docutils literal"><span class="pre">GATGACGGTGT</span></tt> is a 5’
primer sequence we want to look for in some FASTQ formatted read data.
As in the example above, we’ll use the <tt class="docutils literal"><span class="pre">SRR020192.fastq</span></tt> file
downloaded from the ENA
(<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</a>&gt;`__).
The same approach would work with any other supported file format (e.g.
FASTA files).</p>
<p>This code uses <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> with a generator expression (to avoid
loading all the sequences into memory at once), and the <tt class="docutils literal"><span class="pre">Seq</span></tt> object’s
<tt class="docutils literal"><span class="pre">startswith</span></tt> method to see if the read starts with the primer
sequence:</p>
<p>That should find 13819 reads from <tt class="docutils literal"><span class="pre">SRR014849.fastq</span></tt> and save them to a
new FASTQ file, <tt class="docutils literal"><span class="pre">with_primer.fastq</span></tt>.</p>
<p>Now suppose that instead you wanted to make a FASTQ file containing
these reads but with the primer sequence removed? That’s just a small
change as we can slice the <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> (see
Section <a class="reference external" href="#sec:SeqRecord-slicing">4.6</a>) to remove the first eleven
letters (the length of our primer):</p>
<p>Again, that should pull out the 13819 reads from <tt class="docutils literal"><span class="pre">SRR020192.fastq</span></tt>,
but this time strip off the first ten characters, and save them to
another new FASTQ file, <tt class="docutils literal"><span class="pre">with_primer_trimmed.fastq</span></tt>.</p>
<p>Finally, suppose you want to create a new FASTQ file where these reads
have their primer removed, but all the other reads are kept as they
were? If we want to still use a generator expression, it is probably
clearest to define our own trim function:</p>
<p>This takes longer, as this time the output file contains all 41892
reads. Again, we’re used a generator expression to avoid any memory
problems. You could alternatively use a generator function rather than a
generator expression.</p>
<p>This form is more flexible if you want to do something more complicated
where only some of the records are retained – as shown in the next
example.</p>
</div>
<div class="section" id="trimming-off-adaptor-sequences">
<h3>18.1.8  Trimming off adaptor sequences<a class="headerlink" href="#trimming-off-adaptor-sequences" title="Permalink to this headline">¶</a></h3>
<p>This is essentially a simple extension to the previous example. We are
going to going to pretend <tt class="docutils literal"><span class="pre">GATGACGGTGT</span></tt> is an adaptor sequence in some
FASTQ formatted read data, again the <tt class="docutils literal"><span class="pre">SRR020192.fastq</span></tt> file from the
NCBI
(<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</a>&gt;`__).</p>
<p>This time however, we will look for the sequence <em>anywhere</em> in the
reads, not just at the very beginning:</p>
<p>Because we are using a FASTQ input file in this example, the
<tt class="docutils literal"><span class="pre">SeqRecord</span></tt> objects have per-letter-annotation for the quality scores.
By slicing the <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> object the appropriate scores are used on
the trimmed records, so we can output them as a FASTQ file too.</p>
<p>Compared to the output of the previous example where we only looked for
a primer/adaptor at the start of each read, you may find some of the
trimmed reads are quite short after trimming (e.g. if the adaptor was
found in the middle rather than near the start). So, let’s add a minimum
length requirement as well:</p>
<p>By changing the format names, you could apply this to FASTA files
instead. This code also could be extended to do a fuzzy match instead of
an exact match (maybe using a pairwise alignment, or taking into account
the read quality scores), but that will be much slower.</p>
</div>
<div class="section" id="converting-fastq-files">
<h3>18.1.9  Converting FASTQ files<a class="headerlink" href="#converting-fastq-files" title="Permalink to this headline">¶</a></h3>
<p>Back in Section <a class="reference external" href="#sec:SeqIO-conversion">5.5.2</a> we showed how to use
<tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to convert between two file formats. Here we’ll go into a
little more detail regarding FASTQ files which are used in second
generation DNA sequencing. Please refer to Cock <em>et al.</em> (2009)
[<a class="reference external" href="#cock2010">7</a>] for a longer description. FASTQ files store both
the DNA sequence (as a string) and the associated read qualities.</p>
<p>PHRED scores (used in most FASTQ files, and also in QUAL files, ACE
files and SFF files) have become a <em>de facto</em> standard for representing
the probability of a sequencing error (here denoted by <em>P</em><sub>*e*</sub>)
at a given base using a simple base ten log transformation:</p>
<table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><em>Q</em><sub>PHRED</sub> = − 10 × log:sub:<cite>10</cite> ( <em>P</em><sub>*e*</sub> )     (18.1)</td>
</tr>
</tbody>
</table>
<p>This means a wrong read (<em>P</em><sub>*e*</sub> = 1) gets a PHRED quality of 0,
while a very good read like <em>P</em><sub>*e*</sub> = 0.00001 gets a PHRED
quality of 50. While for raw sequencing data qualities higher than this
are rare, with post processing such as read mapping or assembly,
qualities of up to about 90 are possible (indeed, the MAQ tool allows
for PHRED scores in the range 0 to 93 inclusive).</p>
<p>The FASTQ format has the potential to become a <em>de facto</em> standard for
storing the letters and quality scores for a sequencing read in a single
plain text file. The only fly in the ointment is that there are at least
three versions of the FASTQ format which are incompatible and difficult
to distinguish...</p>
<ol class="arabic simple">
<li>The original Sanger FASTQ format uses PHRED qualities encoded with an
ASCII offset of 33. The NCBI are using this format in their Short
Read Archive. We call this the <tt class="docutils literal"><span class="pre">fastq</span></tt> (or <tt class="docutils literal"><span class="pre">fastq-sanger</span></tt>) format
in <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt>.</li>
<li>Solexa (later bought by Illumina) introduced their own version using
Solexa qualities encoded with an ASCII offset of 64. We call this the
<tt class="docutils literal"><span class="pre">fastq-solexa</span></tt> format.</li>
<li>Illumina pipeline 1.3 onwards produces FASTQ files with PHRED
qualities (which is more consistent), but encoded with an ASCII
offset of 64. We call this the <tt class="docutils literal"><span class="pre">fastq-illumina</span></tt> format.</li>
</ol>
<p>The Solexa quality scores are defined using a different log
transformation:</p>
<p><em>Q</em><sub>Solexa</sub> = − 10 × log:sub:<cite>10</cite></p>
<div class="line-block">
<div class="line">⎛</div>
<div class="line-block">
<div class="line">⎜</div>
<div class="line">⎜</div>
<div class="line">⎝</div>
</div>
</div>
<table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><em>P</em><sub>*e*</sub></td>
</tr>
<tr class="row-even"><td>&nbsp;</td>
</tr>
<tr class="row-odd"><td>1−<em>P</em><sub>*e*</sub></td>
</tr>
</tbody>
</table>
<div class="line-block">
<div class="line">⎞</div>
<div class="line-block">
<div class="line">⎟</div>
<div class="line">⎟</div>
<div class="line">⎠</div>
</div>
</div>
<blockquote>
<div>(18.2)</div></blockquote>
<p>Given Solexa/Illumina have now moved to using PHRED scores in version
1.3 of their pipeline, the Solexa quality scores will gradually fall out
of use. If you equate the error estimates (<em>P</em><sub>*e*</sub>) these two
equations allow conversion between the two scoring systems - and
Biopython includes functions to do this in the <tt class="docutils literal"><span class="pre">Bio.SeqIO.QualityIO</span></tt>
module, which are called if you use <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to convert an old
Solexa/Illumina file into a standard Sanger FASTQ file:</p>
<p>If you want to convert a new Illumina 1.3+ FASTQ file, all that gets
changed is the ASCII offset because although encoded differently the
scores are all PHRED qualities:</p>
<p>Note that using <tt class="docutils literal"><span class="pre">Bio.SeqIO.convert()</span></tt> like this is <em>much</em> faster than
combining <tt class="docutils literal"><span class="pre">Bio.SeqIO.parse()</span></tt> and <tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt> because
optimised code is used for converting between FASTQ variants (and also
for FASTQ to FASTA conversion).</p>
<p>For good quality reads, PHRED and Solexa scores are approximately equal,
which means since both the <tt class="docutils literal"><span class="pre">fasta-solexa</span></tt> and <tt class="docutils literal"><span class="pre">fastq-illumina</span></tt>
formats use an ASCII offset of 64 the files are almost the same. This
was a deliberate design choice by Illumina, meaning applications
expecting the old <tt class="docutils literal"><span class="pre">fasta-solexa</span></tt> style files will probably be OK using
the newer <tt class="docutils literal"><span class="pre">fastq-illumina</span></tt> files (on good data). Of course, both
variants are very different from the original FASTQ standard as used by
Sanger, the NCBI, and elsewhere (format name <tt class="docutils literal"><span class="pre">fastq</span></tt> or
<tt class="docutils literal"><span class="pre">fastq-sanger</span></tt>).</p>
<p>For more details, see the built in help (also
<a class="reference external" href="http://www.biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html">online</a>):</p>
</div>
<div class="section" id="converting-fasta-and-qual-files-into-fastq-files">
<h3>18.1.10  Converting FASTA and QUAL files into FASTQ files<a class="headerlink" href="#converting-fasta-and-qual-files-into-fastq-files" title="Permalink to this headline">¶</a></h3>
<p>FASTQ files hold <em>both</em> sequences and their quality strings. FASTA files
hold <em>just</em> sequences, while QUAL files hold <em>just</em> the qualities.
Therefore a single FASTQ file can be converted to or from <em>paired</em> FASTA
and QUAL files.</p>
<p>Going from FASTQ to FASTA is easy:</p>
<p>Going from FASTQ to QUAL is also easy:</p>
<p>However, the reverse is a little more tricky. You can use
<tt class="docutils literal"><span class="pre">Bio.SeqIO.parse()</span></tt> to iterate over the records in a <em>single</em> file,
but in this case we have two input files. There are several strategies
possible, but assuming that the two files are really paired the most
memory efficient way is to loop over both together. The code is a little
fiddly, so we provide a function called <tt class="docutils literal"><span class="pre">PairedFastaQualIterator</span></tt> in
the <tt class="docutils literal"><span class="pre">Bio.SeqIO.QualityIO</span></tt> module to do this. This takes two handles
(the FASTA file and the QUAL file) and returns a <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> iterator:</p>
<p>This function will check that the FASTA and QUAL files are consistent
(e.g. the records are in the same order, and have the same sequence
length). You can combine this with the <tt class="docutils literal"><span class="pre">Bio.SeqIO.write()</span></tt> function to
convert a pair of FASTA and QUAL files into a single FASTQ files:</p>
</div>
<div class="section" id="indexing-a-fastq-file">
<h3>18.1.11  Indexing a FASTQ file<a class="headerlink" href="#indexing-a-fastq-file" title="Permalink to this headline">¶</a></h3>
<p>FASTQ files are often very large, with millions of reads in them. Due to
the sheer amount of data, you can’t load all the records into memory at
once. This is why the examples above (filtering and trimming) iterate
over the file looking at just one <tt class="docutils literal"><span class="pre">SeqRecord</span></tt> at a time.</p>
<p>However, sometimes you can’t use a big loop or an iterator - you may
need random access to the reads. Here the <tt class="docutils literal"><span class="pre">Bio.SeqIO.index()</span></tt> function
may prove very helpful, as it allows you to access any read in the FASTQ
file by its name (see Section <a class="reference external" href="#sec:SeqIO-index">5.4.2</a>).</p>
<p>Again we’ll use the <tt class="docutils literal"><span class="pre">SRR020192.fastq</span></tt> file from the ENA
(<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR020/SRR020192/SRR020192.fastq.gz</a>&gt;`__),
although this is actually quite a small FASTQ file with less than 50,000
reads:</p>
<p>When testing this on a FASTQ file with seven million reads, indexing
took about a minute, but record access was almost instant.</p>
<p>The example in Section <a class="reference external" href="#sec:SeqIO-sort">18.1.5</a> show how you can
use the <tt class="docutils literal"><span class="pre">Bio.SeqIO.index()</span></tt> function to sort a large FASTA file – this
could also be used on FASTQ files.</p>
</div>
<div class="section" id="converting-sff-files">
<h3>18.1.12  Converting SFF files<a class="headerlink" href="#converting-sff-files" title="Permalink to this headline">¶</a></h3>
<p>If you work with 454 (Roche) sequence data, you will probably have
access to the raw data as a Standard Flowgram Format (SFF) file. This
contains the sequence reads (called bases) with quality scores and the
original flow information.</p>
<p>A common task is to convert from SFF to a pair of FASTA and QUAL files,
or to a single FASTQ file. These operations are trivial using the
<tt class="docutils literal"><span class="pre">Bio.SeqIO.convert()</span></tt> function (see
Section <a class="reference external" href="#sec:SeqIO-conversion">5.5.2</a>):</p>
<p>Remember the convert function returns the number of records, in this
example just ten. This will give you the <em>untrimmed</em> reads, where the
leading and trailing poor quality sequence or adaptor will be in lower
case. If you want the <em>trimmed</em> reads (using the clipping information
recorded within the SFF file) use this:</p>
<p>If you run Linux, you could ask Roche for a copy of their “off
instrument” tools (often referred to as the Newbler tools). This offers
an alternative way to do SFF to FASTA or QUAL conversion at the command
line (but currently FASTQ output is not supported), e.g.</p>
<p>The way Biopython uses mixed case sequence strings to represent the
trimming points deliberately mimics what the Roche tools do.</p>
<p>For more information on the Biopython SFF support, consult the built in
help:</p>
</div>
<div class="section" id="identifying-open-reading-frames">
<h3>18.1.13  Identifying open reading frames<a class="headerlink" href="#identifying-open-reading-frames" title="Permalink to this headline">¶</a></h3>
<p>A very simplistic first step at identifying possible genes is to look
for open reading frames (ORFs). By this we mean look in all six frames
for long regions without stop codons – an ORF is just a region of
nucleotides with no in frame stop codons.</p>
<p>Of course, to find a gene you would also need to worry about locating a
start codon, possible promoters – and in Eukaryotes there are introns to
worry about too. However, this approach is still useful in viruses and
Prokaryotes.</p>
<p>To show how you might approach this with Biopython, we’ll need a
sequence to search, and as an example we’ll again use the bacterial
plasmid – although this time we’ll start with a plain FASTA file with no
pre-marked genes:
<tt class="docutils literal"><span class="pre">`NC_005816.fna</span></tt> &lt;<a class="reference external" href="http://biopython.org/SRC/biopython/Tests/GenBank/NC_005816.fna">http://biopython.org/SRC/biopython/Tests/GenBank/NC_005816.fna</a>&gt;`__.
This is a bacterial sequence, so we’ll want to use NCBI codon table 11
(see Section <a class="reference external" href="#sec:translation">3.9</a> about translation).</p>
<p>Here is a neat trick using the <tt class="docutils literal"><span class="pre">Seq</span></tt> object’s <tt class="docutils literal"><span class="pre">split</span></tt> method to get
a list of all the possible ORF translations in the six reading frames:</p>
<p>Note that here we are counting the frames from the 5’ end (start) of
<em>each</em> strand. It is sometimes easier to always count from the 5’ end
(start) of the <em>forward</em> strand.</p>
<p>You could easily edit the above loop based code to build up a list of
the candidate proteins, or convert this to a list comprehension. Now,
one thing this code doesn’t do is keep track of where the proteins are.</p>
<p>You could tackle this in several ways. For example, the following code
tracks the locations in terms of the protein counting, and converts back
to the parent sequence by multiplying by three, then adjusting for the
frame and strand:</p>
<p>And the output:</p>
<p>If you comment out the sort statement, then the protein sequences will
be shown in the same order as before, so you can check this is doing the
same thing. Here we have sorted them by location to make it easier to
compare to the actual annotation in the GenBank file (as visualised in
Section <a class="reference external" href="#sec:gd_nice_example">17.1.9</a>).</p>
<p>If however all you want to find are the locations of the open reading
frames, then it is a waste of time to translate every possible codon,
including doing the reverse complement to search the reverse strand too.
All you need to do is search for the possible stop codons (and their
reverse complements). Using regular expressions is an obvious approach
here (see the Python module <tt class="docutils literal"><span class="pre">re</span></tt>). These are an extremely powerful
(but rather complex) way of describing search strings, which are
supported in lots of programming languages and also command line tools
like <tt class="docutils literal"><span class="pre">grep</span></tt> as well). You can find whole books about this topic!</p>
</div>
</div>
<div class="section" id="sequence-parsing-plus-simple-plots">
<h2>18.2  Sequence parsing plus simple plots<a class="headerlink" href="#sequence-parsing-plus-simple-plots" title="Permalink to this headline">¶</a></h2>
<p>This section shows some more examples of sequence parsing, using the
<tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> module described in Chapter <a class="reference external" href="#chapter:Bio.SeqIO">5</a>,
plus the Python library matplotlib’s <tt class="docutils literal"><span class="pre">pylab</span></tt> plotting interface (see
<a class="reference external" href="http://matplotlib.sourceforge.net/">the matplotlib website for a
tutorial</a>). Note that to follow
these examples you will need matplotlib installed - but without it you
can still try the data parsing bits.</p>
<div class="section" id="histogram-of-sequence-lengths">
<h3>18.2.1  Histogram of sequence lengths<a class="headerlink" href="#histogram-of-sequence-lengths" title="Permalink to this headline">¶</a></h3>
<p>There are lots of times when you might want to visualise the
distribution of sequence lengths in a dataset – for example the range of
contig sizes in a genome assembly project. In this example we’ll reuse
our orchid FASTA file
<tt class="docutils literal"><span class="pre">`ls_orchid.fasta</span></tt> &lt;<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta">http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta</a>&gt;`__
which has only 94 sequences.</p>
<p>First of all, we will use <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to parse the FASTA file and
compile a list of all the sequence lengths. You could do this with a for
loop, but I find a list comprehension more pleasing:</p>
<p>Now that we have the lengths of all the genes (as a list of integers),
we can use the matplotlib histogram function to display it.</p>
<p>That should pop up a new window containing the following graph:</p>
<p><a href="#id1"><span class="problematic" id="id2">|image26|</span></a></p>
<p>Notice that most of these orchid sequences are about 740 bp long, and
there could be two distinct classes of sequence here with a subset of
shorter sequences.</p>
<p><em>Tip:</em> Rather than using <tt class="docutils literal"><span class="pre">pylab.show()</span></tt> to show the plot in a window,
you can also use <tt class="docutils literal"><span class="pre">pylab.savefig(...)</span></tt> to save the figure to a file
(e.g. as a PNG or PDF).</p>
</div>
<div class="section" id="plot-of-sequence-gc">
<h3>18.2.2  Plot of sequence GC%<a class="headerlink" href="#plot-of-sequence-gc" title="Permalink to this headline">¶</a></h3>
<p>Another easily calculated quantity of a nucleotide sequence is the GC%.
You might want to look at the GC% of all the genes in a bacterial genome
for example, and investigate any outliers which could have been recently
acquired by horizontal gene transfer. Again, for this example we’ll
reuse our orchid FASTA file
<tt class="docutils literal"><span class="pre">`ls_orchid.fasta</span></tt> &lt;<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta">http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta</a>&gt;`__.</p>
<p>First of all, we will use <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> to parse the FASTA file and
compile a list of all the GC percentages. Again, you could do this with
a for loop, but I prefer this:</p>
<p>Having read in each sequence and calculated the GC%, we then sorted them
into ascending order. Now we’ll take this list of floating point values
and plot them with matplotlib:</p>
<p>As in the previous example, that should pop up a new window containing a
graph:</p>
<p><a href="#id3"><span class="problematic" id="id4">|image27|</span></a></p>
<p>If you tried this on the full set of genes from one organism, you’d
probably get a much smoother plot than this.</p>
</div>
<div class="section" id="nucleotide-dot-plots">
<h3>18.2.3  Nucleotide dot plots<a class="headerlink" href="#nucleotide-dot-plots" title="Permalink to this headline">¶</a></h3>
<p>A dot plot is a way of visually comparing two nucleotide sequences for
similarity to each other. A sliding window is used to compare short
sub-sequences to each other, often with a mis-match threshold. Here for
simplicity we’ll only look for perfect matches (shown in black in the
plot below).</p>
<p>To start off, we’ll need two sequences. For the sake of argument, we’ll
just take the first two from our orchid FASTA file
<tt class="docutils literal"><span class="pre">`ls_orchid.fasta</span></tt> &lt;<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta">http://biopython.org/DIST/docs/tutorial/examples/ls_orchid.fasta</a>&gt;`__:</p>
<p>We’re going to show two approaches. Firstly, a simple naive
implementation which compares all the window sized sub-sequences to each
other to compiles a similarity matrix. You could construct a matrix or
array object, but here we just use a list of lists of booleans created
with a nested list comprehension:</p>
<p>Note that we have <em>not</em> checked for reverse complement matches here. Now
we’ll use the matplotlib’s <tt class="docutils literal"><span class="pre">pylab.imshow()</span></tt> function to display this
data, first requesting the gray color scheme so this is done in black
and white:</p>
<p>That should pop up a new window containing a graph like this:</p>
<p><a href="#id5"><span class="problematic" id="id6">|image28|</span></a></p>
<p>As you might have expected, these two sequences are very similar with a
partial line of window sized matches along the diagonal. There are no
off diagonal matches which would be indicative of inversions or other
interesting events.</p>
<p>The above code works fine on small examples, but there are two problems
applying this to larger sequences, which we will address below. First
off all, this brute force approach to the all against all comparisons is
very slow. Instead, we’ll compile dictionaries mapping the window sized
sub-sequences to their locations, and then take the set intersection to
find those sub-sequences found in both sequences. This uses more memory,
but is <em>much</em> faster. Secondly, the <tt class="docutils literal"><span class="pre">pylab.imshow()</span></tt> function is
limited in the size of matrix it can display. As an alternative, we’ll
use the <tt class="docutils literal"><span class="pre">pylab.scatter()</span></tt> function.</p>
<p>We start by creating dictionaries mapping the window-sized sub-sequences
to locations:</p>
<p>In order to use the <tt class="docutils literal"><span class="pre">pylab.scatter()</span></tt> we need separate lists for the
<em>x</em> and <em>y</em> co-ordinates:</p>
<p>We are now ready to draw the revised dot plot as a scatter plot:</p>
<p>That should pop up a new window containing a graph like this:</p>
<p><a href="#id7"><span class="problematic" id="id8">|image29|</span></a></p>
<p>Personally I find this second plot much easier to read! Again note that
we have <em>not</em> checked for reverse complement matches here – you could
extend this example to do this, and perhaps plot the forward matches in
one color and the reverse matches in another.</p>
</div>
<div class="section" id="plotting-the-quality-scores-of-sequencing-read-data">
<h3>18.2.4  Plotting the quality scores of sequencing read data<a class="headerlink" href="#plotting-the-quality-scores-of-sequencing-read-data" title="Permalink to this headline">¶</a></h3>
<p>If you are working with second generation sequencing data, you may want
to try plotting the quality data. Here is an example using two FASTQ
files containing paired end reads, <tt class="docutils literal"><span class="pre">SRR001666_1.fastq</span></tt> for the forward
reads, and <tt class="docutils literal"><span class="pre">SRR001666_2.fastq</span></tt> for the reverse reads. These were
downloaded from the ENA sequence read archive FTP site
(<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_1.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_1.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_1.fastq.gz</a>&gt;`__
and
<tt class="docutils literal"><span class="pre">`ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_2.fastq.gz</span></tt> &lt;<a class="reference external" href="ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_2.fastq.gz">ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR001/SRR001666/SRR001666_2.fastq.gz</a>&gt;`__),
and are from <em>E. coli</em> – see
<tt class="docutils literal"><span class="pre">`http://www.ebi.ac.uk/ena/data/view/SRR001666</span></tt> &lt;<a class="reference external" href="http://www.ebi.ac.uk/ena/data/view/SRR001666">http://www.ebi.ac.uk/ena/data/view/SRR001666</a>&gt;`__
for details. In the following code the <tt class="docutils literal"><span class="pre">pylab.subplot(...)</span></tt> function
is used in order to show the forward and reverse qualities on two
subplots, side by side. There is also a little bit of code to only plot
the first fifty reads.</p>
<p>You should note that we are using the <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt> format name
<tt class="docutils literal"><span class="pre">fastq</span></tt> here because the NCBI has saved these reads using the standard
Sanger FASTQ format with PHRED scores. However, as you might guess from
the read lengths, this data was from an Illumina Genome Analyzer and was
probably originally in one of the two Solexa/Illumina FASTQ variant file
formats instead.</p>
<p>This example uses the <tt class="docutils literal"><span class="pre">pylab.savefig(...)</span></tt> function instead of
<tt class="docutils literal"><span class="pre">pylab.show(...)</span></tt>, but as mentioned before both are useful. Here is
the result:</p>
<p><a href="#id9"><span class="problematic" id="id10">|image30|</span></a></p>
</div>
</div>
<div class="section" id="dealing-with-alignments">
<h2>18.3  Dealing with alignments<a class="headerlink" href="#dealing-with-alignments" title="Permalink to this headline">¶</a></h2>
<p>This section can been seen as a follow on to
Chapter <a class="reference external" href="#chapter:Bio.AlignIO">6</a>.</p>
<div class="section" id="calculating-summary-information">
<h3>18.3.1  Calculating summary information<a class="headerlink" href="#calculating-summary-information" title="Permalink to this headline">¶</a></h3>
<p>Once you have an alignment, you are very likely going to want to find
out information about it. Instead of trying to have all of the functions
that can generate information about an alignment in the alignment object
itself, we’ve tried to separate out the functionality into separate
classes, which act on the alignment.</p>
<p>Getting ready to calculate summary information about an object is quick
to do. Let’s say we’ve got an alignment object called <tt class="docutils literal"><span class="pre">alignment</span></tt>, for
example read in using <tt class="docutils literal"><span class="pre">Bio.AlignIO.read(...)</span></tt> as described in
Chapter <a class="reference external" href="#chapter:Bio.AlignIO">6</a>. All we need to do to get an
object that will calculate summary information is:</p>
<p>The <tt class="docutils literal"><span class="pre">summary_align</span></tt> object is very useful, and will do the following
neat things for you:</p>
<ol class="arabic simple">
<li>Calculate a quick consensus sequence – see
section <a class="reference external" href="#sec:consensus">18.3.2</a></li>
<li>Get a position specific score matrix for the alignment – see
section <a class="reference external" href="#sec:pssm">18.3.3</a></li>
<li>Calculate the information content for the alignment – see
section <a class="reference external" href="#sec:getting_info_content">18.3.4</a></li>
<li>Generate information on substitutions in the alignment –
section <a class="reference external" href="#sec:sub_matrix">18.4</a> details using this to generate a
substitution matrix.</li>
</ol>
</div>
<div class="section" id="calculating-a-quick-consensus-sequence">
<h3>18.3.2  Calculating a quick consensus sequence<a class="headerlink" href="#calculating-a-quick-consensus-sequence" title="Permalink to this headline">¶</a></h3>
<p>The <tt class="docutils literal"><span class="pre">SummaryInfo</span></tt> object, described in
section <a class="reference external" href="#sec:summary_info">18.3.1</a>, provides functionality to
calculate a quick consensus of an alignment. Assuming we’ve got a
<tt class="docutils literal"><span class="pre">SummaryInfo</span></tt> object called <tt class="docutils literal"><span class="pre">summary_align</span></tt> we can calculate a
consensus by doing:</p>
<p>As the name suggests, this is a really simple consensus calculator, and
will just add up all of the residues at each point in the consensus, and
if the most common value is higher than some threshold value will add
the common residue to the consensus. If it doesn’t reach the threshold,
it adds an ambiguity character to the consensus. The returned consensus
object is Seq object whose alphabet is inferred from the alphabets of
the sequences making up the consensus. So doing a <tt class="docutils literal"><span class="pre">print</span> <span class="pre">consensus</span></tt>
would give:</p>
<p>You can adjust how <tt class="docutils literal"><span class="pre">dumb_consensus</span></tt> works by passing optional
parameters:</p>
<blockquote>
<div><dl class="docutils">
<dt><strong>the threshold</strong></dt>
<dd>This is the threshold specifying how common a particular residue has
to be at a position before it is added. The default is 0.7 (meaning
70%).</dd>
</dl>
</div></blockquote>
<dl class="docutils">
<dt><strong>the ambiguous character</strong></dt>
<dd>This is the ambiguity character to use. The default is ’N’.</dd>
<dt><strong>the consensus alphabet</strong></dt>
<dd>This is the alphabet to use for the consensus sequence. If an
alphabet is not specified than we will try to guess the alphabet
based on the alphabets of the sequences in the alignment.</dd>
</dl>
</div>
<div class="section" id="position-specific-score-matrices">
<h3>18.3.3  Position Specific Score Matrices<a class="headerlink" href="#position-specific-score-matrices" title="Permalink to this headline">¶</a></h3>
<p>Position specific score matrices (PSSMs) summarize the alignment
information in a different way than a consensus, and may be useful for
different tasks. Basically, a PSSM is a count matrix. For each column in
the alignment, the number of each alphabet letters is counted and
totaled. The totals are displayed relative to some representative
sequence along the left axis. This sequence may be the consesus
sequence, but can also be any sequence in the alignment. For instance
for the alignment,</p>
<p>the PSSM is:</p>
<p>Let’s assume we’ve got an alignment object called <tt class="docutils literal"><span class="pre">c_align</span></tt>. To get a
PSSM with the consensus sequence along the side we first get a summary
object and calculate the consensus sequence:</p>
<p>Now, we want to make the PSSM, but ignore any <tt class="docutils literal"><span class="pre">N</span></tt> ambiguity residues
when calculating this:</p>
<p>Two notes should be made about this:</p>
<ol class="arabic simple">
<li>To maintain strictness with the alphabets, you can only include
characters along the top of the PSSM that are in the alphabet of the
alignment object. Gaps are not included along the top axis of the
PSSM.</li>
<li>The sequence passed to be displayed along the left side of the axis
does not need to be the consensus. For instance, if you wanted to
display the second sequence in the alignment along this axis, you
would need to do:</li>
</ol>
<p>The command above returns a <tt class="docutils literal"><span class="pre">PSSM</span></tt> object. To print out the PSSM as we
showed above, we simply need to do a <tt class="docutils literal"><span class="pre">print</span> <span class="pre">my_pssm</span></tt>, which gives:</p>
<p>You can access any element of the PSSM by subscripting like
<tt class="docutils literal"><span class="pre">your_pssm[sequence_number][residue_count_name]</span></tt>. For instance, to get
the counts for the ’A’ residue in the second element of the above PSSM
you would do:</p>
<p>The structure of the PSSM class hopefully makes it easy both to access
elements and to pretty print the matrix.</p>
</div>
<div class="section" id="information-content">
<h3>18.3.4  Information Content<a class="headerlink" href="#information-content" title="Permalink to this headline">¶</a></h3>
<p>A potentially useful measure of evolutionary conservation is the
information content of a sequence.</p>
<p>A useful introduction to information theory targeted towards molecular
biologists can be found at
<tt class="docutils literal"><span class="pre">`http://www.lecb.ncifcrf.gov/~toms/paper/primer/</span></tt> &lt;<a class="reference external" href="http://www.lecb.ncifcrf.gov/~toms/paper/primer/">http://www.lecb.ncifcrf.gov/~toms/paper/primer/</a>&gt;`__.
For our purposes, we will be looking at the information content of a
consesus sequence, or a portion of a consensus sequence. We calculate
information content at a particular column in a multiple sequence
alignment using the following formula:</p>
<p><em>IC</em><sub>*j*</sub> =</p>
<table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><em>N</em><sub>*a*</sub></td>
</tr>
<tr class="row-even"><td>∑</td>
</tr>
<tr class="row-odd"><td><em>i</em>=1</td>
</tr>
</tbody>
</table>
<blockquote>
<div><em>P</em><sub>*ij*</sub> <em>log</em></div></blockquote>
<div class="line-block">
<div class="line">⎛</div>
<div class="line-block">
<div class="line">⎜</div>
<div class="line">⎜</div>
<div class="line">⎝</div>
</div>
</div>
<table border="1" class="docutils">
<colgroup>
<col width="100%" />
</colgroup>
<tbody valign="top">
<tr class="row-odd"><td><em>P</em><sub>*ij*</sub></td>
</tr>
<tr class="row-even"><td>&nbsp;</td>
</tr>
<tr class="row-odd"><td><em>Q</em><sub>*i*</sub></td>
</tr>
</tbody>
</table>
<div class="line-block">
<div class="line">⎞</div>
<div class="line-block">
<div class="line">⎟</div>
<div class="line">⎟</div>
<div class="line">⎠</div>
</div>
</div>
<p>where:</p>
<ul class="simple">
<li><em>IC</em><sub>*j*</sub> – The information content for the <em>j</em>-th column in
an alignment.</li>
<li><em>N</em><sub>*a*</sub> – The number of letters in the alphabet.</li>
<li><em>P</em><sub>*ij*</sub> – The frequency of a particular letter <em>i</em> in the
<em>j</em>-th column (i. e. if G occurred 3 out of 6 times in an aligment
column, this would be 0.5)</li>
<li><em>Q</em><sub>*i*</sub> – The expected frequency of a letter <em>i</em>. This is an
optional argument, usage of which is left at the user’s discretion.
By default, it is automatically assigned to 0.05 = 1/20 for a protein
alphabet, and 0.25 = 1/4 for a nucleic acid alphabet. This is for
geting the information content without any assumption of prior
distributions. When assuming priors, or when using a non-standard
alphabet, you should supply the values for <em>Q</em><sub>*i*</sub>.</li>
</ul>
<p>Well, now that we have an idea what information content is being
calculated in Biopython, let’s look at how to get it for a particular
region of the alignment.</p>
<p>First, we need to use our alignment to get an alignment summary object,
which we’ll assume is called <tt class="docutils literal"><span class="pre">summary_align</span></tt> (see
section <a class="reference external" href="#sec:summary_info">18.3.1</a>) for instructions on how to get
this. Once we’ve got this object, calculating the information content
for a region is as easy as:</p>
<p>Wow, that was much easier then the formula above made it look! The
variable <tt class="docutils literal"><span class="pre">info_content</span></tt> now contains a float value specifying the
information content over the specified region (from 5 to 30 of the
alignment). We specifically ignore the ambiguity residue ’N’ when
calculating the information content, since this value is not included in
our alphabet (so we shouldn’t be interested in looking at it!).</p>
<p>As mentioned above, we can also calculate relative information content
by supplying the expected frequencies:</p>
<p>The expected should not be passed as a raw dictionary, but instead by
passed as a <tt class="docutils literal"><span class="pre">SubsMat.FreqTable</span></tt> object (see
section <a class="reference external" href="#sec:freq_table">20.2.2</a> for more information about
FreqTables). The FreqTable object provides a standard for associating
the dictionary with an Alphabet, similar to how the Biopython Seq class
works.</p>
<p>To create a FreqTable object, from the frequency dictionary you just
need to do:</p>
<p>Now that we’ve got that, calculating the relative information content
for our region of the alignment is as simple as:</p>
<p>Now, <tt class="docutils literal"><span class="pre">info_content</span></tt> will contain the relative information content over
the region in relation to the expected frequencies.</p>
<p>The value return is calculated using base 2 as the logarithm base in the
formula above. You can modify this by passing the parameter <tt class="docutils literal"><span class="pre">log_base</span></tt>
as the base you want:</p>
<p>Well, now you are ready to calculate information content. If you want to
try applying this to some real life problems, it would probably be best
to dig into the literature on information content to get an idea of how
it is used. Hopefully your digging won’t reveal any mistakes made in
coding this function!</p>
</div>
</div>
<div class="section" id="substitution-matrices">
<h2>18.4  Substitution Matrices<a class="headerlink" href="#substitution-matrices" title="Permalink to this headline">¶</a></h2>
<p>Substitution matrices are an extremely important part of everyday
bioinformatics work. They provide the scoring terms for classifying how
likely two different residues are to substitute for each other. This is
essential in doing sequence comparisons. The book “Biological Sequence
Analysis” by Durbin et al. provides a really nice introduction to
Substitution Matrices and their uses. Some famous substitution matrices
are the PAM and BLOSUM series of matrices.</p>
<p>Biopython provides a ton of common substitution matrices, and also
provides functionality for creating your own substitution matrices.</p>
<div class="section" id="using-common-substitution-matrices">
<h3>18.4.1  Using common substitution matrices<a class="headerlink" href="#using-common-substitution-matrices" title="Permalink to this headline">¶</a></h3>
</div>
<div class="section" id="creating-your-own-substitution-matrix-from-an-alignment">
<h3>18.4.2  Creating your own substitution matrix from an alignment<a class="headerlink" href="#creating-your-own-substitution-matrix-from-an-alignment" title="Permalink to this headline">¶</a></h3>
<p>A very cool thing that you can do easily with the substitution matrix
classes is to create your own substitution matrix from an alignment. In
practice, this is normally done with protein alignments. In this
example, we’ll first get a Biopython alignment object and then get a
summary object to calculate info about the alignment. The file
containing <a class="reference external" href="examples/protein.aln">protein.aln</a> (also available online
<a class="reference external" href="http://biopython.org/DIST/docs/tutorial/examples/protein.aln">here</a>)
contains the Clustalw alignment output.</p>
<p>Sections <a class="reference external" href="#sec:align_clustal">6.4.1</a>
and <a class="reference external" href="#sec:summary_info">18.3.1</a> contain more information on doing
this.</p>
<p>Now that we’ve got our <tt class="docutils literal"><span class="pre">summary_align</span></tt> object, we want to use it to
find out the number of times different residues substitute for each
other. To make the example more readable, we’ll focus on only amino
acids with polar charged side chains. Luckily, this can be done easily
when generating a replacement dictionary, by passing in all of the
characters that should be ignored. Thus we’ll create a dictionary of
replacements for only charged polar amino acids using:</p>
<p>This information about amino acid replacements is represented as a
python dictionary which will look something like (the order can vary):</p>
<p>This information gives us our accepted number of replacements, or how
often we expect different things to substitute for each other. It turns
out, amazingly enough, that this is all of the information we need to go
ahead and create a substitution matrix. First, we use the replacement
dictionary information to create an Accepted Replacement Matrix (ARM):</p>
<p>With this accepted replacement matrix, we can go right ahead and create
our log odds matrix (i. e. a standard type Substitution Matrix):</p>
<p>The log odds matrix you create is customizable with the following
optional arguments:</p>
<ul class="simple">
<li><tt class="docutils literal"><span class="pre">exp_freq_table</span></tt> – You can pass a table of expected frequencies for
each alphabet. If supplied, this will be used instead of the passed
accepted replacement matrix when calculate expected replacments.</li>
<li><tt class="docutils literal"><span class="pre">logbase</span></tt> - The base of the logarithm taken to create the log odd
matrix. Defaults to base 10.</li>
<li><tt class="docutils literal"><span class="pre">factor</span></tt> - The factor to multiply each matrix entry by. This
defaults to 10, which normally makes the matrix numbers easy to work
with.</li>
<li><tt class="docutils literal"><span class="pre">round_digit</span></tt> - The digit to round to in the matrix. This defaults
to 0 (i. e. no digits).</li>
</ul>
<p>Once you’ve got your log odds matrix, you can display it prettily using
the function <tt class="docutils literal"><span class="pre">print_mat</span></tt>. Doing this on our created matrix gives:</p>
<p>Very nice. Now we’ve got our very own substitution matrix to play with!</p>
</div>
</div>
<div class="section" id="biosql-storing-sequences-in-a-relational-database">
<h2>18.5  BioSQL – storing sequences in a relational database<a class="headerlink" href="#biosql-storing-sequences-in-a-relational-database" title="Permalink to this headline">¶</a></h2>
<p><a class="reference external" href="http://www.biosql.org/">BioSQL</a> is a joint effort between the
<a class="reference external" href="http://open-bio.org/">OBF</a> projects (BioPerl, BioJava etc) to
support a shared database schema for storing sequence data. In theory,
you could load a GenBank file into the database with BioPerl, then using
Biopython extract this from the database as a record object with
features - and get more or less the same thing as if you had loaded the
GenBank file directly as a SeqRecord using <tt class="docutils literal"><span class="pre">Bio.SeqIO</span></tt>
(Chapter <a class="reference external" href="#chapter:Bio.SeqIO">5</a>).</p>
<p>Biopython’s BioSQL module is currently documented at
<tt class="docutils literal"><span class="pre">`http://biopython.org/wiki/BioSQL</span></tt> &lt;<a class="reference external" href="http://biopython.org/wiki/BioSQL">http://biopython.org/wiki/BioSQL</a>&gt;`__
which is part of our wiki pages.</p>
</div>
</div>


          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
  <h3><a href="index.html">Table Of Contents</a></h3>
  <ul>
<li><a class="reference internal" href="#">Chapter 18  Cookbook – Cool things to do with it</a><ul>
<li><a class="reference internal" href="#working-with-sequence-files">18.1  Working with sequence files</a><ul>
<li><a class="reference internal" href="#filtering-a-sequence-file">18.1.1  Filtering a sequence file</a></li>
<li><a class="reference internal" href="#producing-randomised-genomes">18.1.2  Producing randomised genomes</a></li>
<li><a class="reference internal" href="#translating-a-fasta-file-of-cds-entries">18.1.3  Translating a FASTA file of CDS entries</a></li>
<li><a class="reference internal" href="#making-the-sequences-in-a-fasta-file-upper-case">18.1.4  Making the sequences in a FASTA file upper case</a></li>
<li><a class="reference internal" href="#sorting-a-sequence-file">18.1.5  Sorting a sequence file</a></li>
<li><a class="reference internal" href="#simple-quality-filtering-for-fastq-files">18.1.6  Simple quality filtering for FASTQ files</a></li>
<li><a class="reference internal" href="#trimming-off-primer-sequences">18.1.7  Trimming off primer sequences</a></li>
<li><a class="reference internal" href="#trimming-off-adaptor-sequences">18.1.8  Trimming off adaptor sequences</a></li>
<li><a class="reference internal" href="#converting-fastq-files">18.1.9  Converting FASTQ files</a></li>
<li><a class="reference internal" href="#converting-fasta-and-qual-files-into-fastq-files">18.1.10  Converting FASTA and QUAL files into FASTQ files</a></li>
<li><a class="reference internal" href="#indexing-a-fastq-file">18.1.11  Indexing a FASTQ file</a></li>
<li><a class="reference internal" href="#converting-sff-files">18.1.12  Converting SFF files</a></li>
<li><a class="reference internal" href="#identifying-open-reading-frames">18.1.13  Identifying open reading frames</a></li>
</ul>
</li>
<li><a class="reference internal" href="#sequence-parsing-plus-simple-plots">18.2  Sequence parsing plus simple plots</a><ul>
<li><a class="reference internal" href="#histogram-of-sequence-lengths">18.2.1  Histogram of sequence lengths</a></li>
<li><a class="reference internal" href="#plot-of-sequence-gc">18.2.2  Plot of sequence GC%</a></li>
<li><a class="reference internal" href="#nucleotide-dot-plots">18.2.3  Nucleotide dot plots</a></li>
<li><a class="reference internal" href="#plotting-the-quality-scores-of-sequencing-read-data">18.2.4  Plotting the quality scores of sequencing read data</a></li>
</ul>
</li>
<li><a class="reference internal" href="#dealing-with-alignments">18.3  Dealing with alignments</a><ul>
<li><a class="reference internal" href="#calculating-summary-information">18.3.1  Calculating summary information</a></li>
<li><a class="reference internal" href="#calculating-a-quick-consensus-sequence">18.3.2  Calculating a quick consensus sequence</a></li>
<li><a class="reference internal" href="#position-specific-score-matrices">18.3.3  Position Specific Score Matrices</a></li>
<li><a class="reference internal" href="#information-content">18.3.4  Information Content</a></li>
</ul>
</li>
<li><a class="reference internal" href="#substitution-matrices">18.4  Substitution Matrices</a><ul>
<li><a class="reference internal" href="#using-common-substitution-matrices">18.4.1  Using common substitution matrices</a></li>
<li><a class="reference internal" href="#creating-your-own-substitution-matrix-from-an-alignment">18.4.2  Creating your own substitution matrix from an alignment</a></li>
</ul>
</li>
<li><a class="reference internal" href="#biosql-storing-sequences-in-a-relational-database">18.5  BioSQL – storing sequences in a relational database</a></li>
</ul>
</li>
</ul>

  <h4>Previous topic</h4>
  <p class="topless"><a href="chr17.html"
                        title="previous chapter">Chapter 17  Graphics including GenomeDiagram</a></p>
  <h4>Next topic</h4>
  <p class="topless"><a href="chr19.html"
                        title="next chapter">Chapter 19  The Biopython testing framework</a></p>
  <h3>This Page</h3>
  <ul class="this-page-menu">
    <li><a href="_sources/chr18.txt"
           rel="nofollow">Show Source</a></li>
  </ul>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="chr19.html" title="Chapter 19 The Biopython testing framework"
             >next</a> |</li>
        <li class="right" >
          <a href="chr17.html" title="Chapter 17 Graphics including GenomeDiagram"
             >previous</a> |</li>
        <li><a href="index.html">Biopython_en 1.0 documentation</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2013, Biopython.
      Created using <a href="http://sphinx-doc.org/">Sphinx</a> 1.2b1.
    </div>
  </body>
</html>